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LESSONS  LEARNED  FROM  RANS  SIMULATIONS 
OF  SHOCK- WAVE/BOUND  ARY-LAYER  INTERACTIONS 
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Usual  two-equation  turbulence  models,  based  on  a  constant  value  of  in  the  de¬ 
finition  of  the  eddy  viscosity,  are  known  to  fail  in  predicting  most  features  of  Shock- 
Wave /Boundary-Layer  Interactions  (SWBLI).  In  searching  a  simple  way  to  improve  their 
capabilities,  it  appears  that  three  different  theoretical  approaches  end  up  with  simple 
and  similar  Weakly  Non-Linear  (WNL)  corrections.  Such  a  correction,  originally  pro¬ 
posed  to  ensure  realizability  in  subsonic  flows,  is  extended  to  deal  with  compressible 
flows  and  applied  to  three-dimensional  supersonic  SWBLI  developing  on  three  single 
sharp  fin  plate  configurations.  The  numerical  solutions  are  obtained  by  solving  the  full 
Reynolds- Averaged  Navier-Stokes  equations  on  grids  up  to  3.3  million  cells  with  the  lin¬ 
ear  and  WNL  versions  of  the  k  —  U)  turbulence  model.  The  WNL  correction  allows  full 
grid-convergence  and  yields  much  better  numerical  solutions,  including  the  prediction 
of  the  pressure  plateau  under  the  lambda  foot  of  the  shock,  the  maximum  skin-friction 
coefficient  on  the  bottom  plate,  and  the  secondary  separation  in  the  strong  SWBLI  case. 
The  improvement  is  due  to  halving  the  turbulence  intensity  in  the  vortical  flow  embedded 
within  the  lambda  foot  of  the  shock.  The  WNL  correction  has  also  a  beneficial  effect  in 
transonic  cases,  as  illustrated  for  a  two-dimensional  channel  flow  over  a  bump.  Finally,  its 
success  in  predicting  unsteady  features  in  the  case  of  shock-induced  oscillations  in  tran¬ 
sonic  flows  over  airfoils  is  underlined.  In  the  considered  3-D  supersonic,  2-D  transonic 
steady  and  unsteady  SWBLI,  the  key-point  is  the  dependency  of  on  the  strain  and 
vorticity  invariants  rather  than  the  non-linear  expansion  of  the  shear  stress. 


Introduction 

One  of  the  main  issues  faced  by  turbulence  mod¬ 
eling  is  to  predict  the  behavior  of  flows  in  presence 
of  large  pressure  gradients,  in  particular  when  yield¬ 
ing  a  boundary-layer  separation.  Because  they  ap¬ 
pear  in  many  aeronautical  applications,  and  because 
they  include  very  large  pressure  gradients,  Shock- 
Wave/Boundary-Layer  Interactions  (SWBLI)  are  a 
very  interesting  archetype  of  such  flows.  In  a  com¬ 
prehensive  and  recent  review,  Dolling1  underlines  the 
main  features  of  these  flows  that  cannot  be  predicted 
today:  heat  transfer,  skin  friction  and  unsteady  pres¬ 
sure  loads.  The  most  recent  monographs  on  the  sub¬ 
ject  are  due  to  Delery  &  Panaras2  and  Smits  &  Dus- 
sauge.3  The  last  review  of  computational  capabilities 
for  SWBLI  is  due  to  Knight  &  Degrez,4  and  the  next 
one  will  come  from  the  RTO/AVT  Working  Group  10. 
Specifically  for  sharp  fin  plate  interactions,  recent  and 
thorough  analyses  can  be  found  in  Refs.  5-7  for  single 
fins,  and  Refs.  8-13  for  double  fins. 

Usual  turbulence  models,  among  which  most  of  two- 
equation  turbulence  models,  are  based  on  the  lin¬ 
ear  formulation  of  the  Boussinesq  assumption:  the 
Reynolds  stress  is  directly  proportional  to  the  mean 
strain  rate.  In  the  following,  such  models  are  called  lin¬ 
ear,  while  models  not  observing  the  linear  Boussinesq 
assumption  are  called  non-linear.  Linear  models  usu¬ 
ally  fail  in  predicting  the  boundary-layer  separation 
or  the  level  of  the  turbulent  kinetic  energy  (and  con- 
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sequently  skin  friction  and  heat  transfer)  in  SWBLI. 
To  overcome  the  difficulties  encountered  by  Reynolds- 
averaged  Navier-Stokes  (RANS)  solvers,  Knight  &  De¬ 
grez4  recommend  developing  Large-Eddy  Simulation 
solvers.  This  is  currently  carried  out,  with  appli¬ 
cations  to  two-dimensional  compressible  high-speed 
flows. 14-16  However,  the  way  still  seems  to  be  long 
before  these  techniques  become  reliable  and  affordable 
in  complex  three-dimensional  configurations.  On  the 
RANS  side,  a  lot  of  efforts  is  put  in  deriving  non¬ 
linear  turbulence  models,  either  in  an  explicit  algebraic 
form  or  through  transport  equations  for  the  Reynolds 
stress  components.  These  models  contain  many  em¬ 
pirical  constants  and  damping  functions  near  walls, 
which  are  not  easy  to  determine  individually  from  a 
set  of  experiments.  This  on-going  research  has  not 
yet  reached  a  point  where  a  unique  formulation  could 
be  able  to  model  a  very  large  range  of  flow  configu¬ 
rations.  In  consequence,  Non-Linear  Eddy- Viscosity 
Models  or  Reynolds-Stress  Transport-Equation  Mod¬ 
els  still  belong  to  the  research  community  and  are 
not  yet  used  intensively  for  engineering  applications. 
The  aim  of  this  paper  is  to  point  out  some  intermedi¬ 
ate  model,  between  usual  linear  two-equation  models 
and  fully  non-linear  models,  which  would  yield  better 
predictions  of  SWBLI  than  linear  models,  and  would 
be  much  simpler  and  more  reliable  than  current  non¬ 
linear  models.  The  idea  is  to  identify  a  correction 
to  linear  models  which  could  be  viewed  as  the  first- 
order  term  of  their  difference  to  non-linear  models. 
Such  a  correction  has  been  recently  pointed  out  af¬ 
ter  a  parametric  study  of  realizability  corrections.17 
It  is  proven  to  be  useful  to  improve  the  prediction  of 
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wall  heat  fluxes  in  Crossing  SWBLI.  Moreover,  it  is 
shown  that  no  other  usual  compressibility  correction 
or  turbulent  lenght-scale  limit  can  further  improve  the 
results.  Here,  different  ways  of  building  Weakly  Non- 
Linear  (WNL)  corrections  are  examined,  together  with 
the  applications  they  have  been  designed  for,  in  order 
to  underline  the  currently  known  applicability  of  the 
proposed  correction.  Then  a  new  set  of  computations 
about  single  sharp  fin  plate  SWBLI  is  analyzed  to  as¬ 
sess  the  ability  of  the  WNL  correction  to  solve  the 
SWBLI  skin-friction  issue.  Finally,  the  effect  of  the 
correction  is  presented  for  steady  and  unsteady  tran¬ 
sonic  interactions. 


Ways  Toward  Weakly  Non-Linear 
Turbulence  Models 

Linear  Models  and  the  Bradshaw’s  Assumption 

Usual  turbulence  models  are  based  on  a  linear  rela¬ 
tionship  between  the  turbulent  stresses  and  the  mean 
strain  rate  derived  from  the  Boussinesq’s  assumption. 
If  ul  is  the  fluctuation  of  the  mass-averaged  velocity 
field  Ut,  the  deviation  of  the  turbulent  stress  writes: 


pub'll j  T  ^  pk  —  2  pf 


(1) 


where  the  mean  strain  is: 

2  \dx0  +  dxj 


(2) 


The  eddy  viscosity  Ut  depends  only  on  turbulent 
scales,  such  as  the  turbulent  kinetic  energy  k  (TKE) 
and  its  specific  dissipation  c o: 


k 

k't  —  P  Cf_i 

U) 


(3) 


Usually,  is  set  to  the  constant  value  c°  =  0.09  which 
is  adequate  to  model  the  flat  plate  boundary  layer. 
However,  this  assumption  is  known  to  be  highly  unre¬ 
alistic  in  presence  of  adverse  pressure  gradients.  In  a 
two-dimensional  boundary  layer,  the  shear  stress  is  no 
more  proportional  to  the  velocity  gradient  but  to  the 
TKE: 

-uv  =  y/<%  k  (4) 

as  pointed  out  by  Bradshaw.18  The  simplest  frame- 
independent  way  to  reconcile  both  expressions  is  to 
introduce  a  weak  non-linearity  in  the  eddy  viscosity, 
such  as: 


which  involves  the  dimensionless  mean  strain  rate: 

s  =  -  with  S2  =  2 Sl3S3%  -  \s2kk  (6) 

CO  o 


This  idea,  originally  proposed  by  Coakley  for  incom¬ 
pressible  flows,19  was  used  again  by  Menter20  to  define 


the  k  —  ui  Shear-Stress  Transport  (SST)  model.  The 
only  difference  in  is  that,  in  extending  the  velocity 
gradient  dU/dy  to  three  dimensions,  Menter  chose  the 
vorticity  invariant: 


2  i  ,  j  i  1  j 


10 


where  fh,  =  „ 
2 


1  fdU,  dU , 


dx ,  dx. 


in  Eq.  (5)  instead  of  the  strain  rate  invariant  s. 


(7) 


Link  with  Realizability 

A  well-known  example  of  problems  encountered  in 
flows  with  strong  velocity  gradients  (and  streamline 
curvature)  is  the  stagnation-point  anomaly:  near  the 
stagnation  point  of  an  impinging  jet21  or  a  turbine 
blade,22  linear  two-equation  models  overstimate  the 
TKE  level  by  orders  of  magnitude.  In  the  case  of 
an  impinging  jet,  the  heat-transfer  coefficient  is  sys¬ 
tematically  overestimated.  Based  on  the  realizability 
principle,23  a  minimal  correction  was  derived  for  two- 
equation  turbulence  models  and  was  shown  to  cure 
the  stagnation-point  anomaly21  and  to  solve  the  jet- 
impingement  problem.24  Although  starting  from  the 
theoretical  realizability  principle  rather  than  the  em¬ 
pirical  Bradshaw’s  assumption,  Durbin  ends  up  with 
a  correction  which  is  formally  identical  to  Coakley’s 
proposal: 


=  min  jc°,  -^=j  with  c= ^  (8) 

for  which  the  constant  cj\J 3  =  0.29  is  very  close  to 
y/dji  =  0.30.*  More  elaborate  expressions  of  have 
been  derived  to  respect  the  realizability  principle,  in¬ 
cluding  formula  where  \/(s2  +  w2)/2  is  substituted 
to  s  in  Eq.  (8)  in  order  to  account  for  some  three- 
dimensional  and  rotational  effects.22  In  most  cases, 
including  those  considered  here,  both  expressions  are 
equivalent.  Whatever  the  expression,  the  essential 
point  is  to  reproduce  the  asymptotic  behavior  of 
when  s  or  \/(s2  -f  w2)/2  tends  toward  infinity. 

The  model  is  completed  with  the  transport  equa¬ 
tions: 


A  (pk)  =  Pk 


frpk0J+^3 


(p  +  a*pt) 


dk 

dx,  _ 


(9) 


A  (pio)  =  a 


k  \  a 


+  Phc  )  —  PpOJ  + 


A 


dx i 


(p,  +  apt ) 


dio 

dxq 


(10) 

where  the  constants  from  the  original  model  are:  a  = 
5/9,  ji  =  5/6,  /3*  =  1,  a  =  a*  =  1/2,  and  the  TKE 
production  rate  is: 


Pk  =  ptS2  ~^Pk  Skk  (11) 

_  Phi  Phc 

"In  Eq.  (8),  only  c  <  1  is  needed  to  ensure  realizability;  the 
value  c  =  1/2  is  required  to  obtain  a  better  fit  of  heat-transfer 
coefficient  measurements  in  impinging  jets.24 
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Different  ways  may  be  found  to  write  the  net  produc¬ 
tion  rate  of  the  dissipation.  Coakley19  and  Moore 
&  Moore22  do  not  modify  the  formal  expression  of 
the  transport  equation  for  the  dissipation  (a„  =  1  in 
Eq.  (10)).  Consequently,  the  decrease  in  Ht  propa¬ 
gates  to  both  TKE  and  dissipation  production  rates 
through  Eq.  (11).  Menter  in  the  SST  model20  and 
Durbin21  divide  the  production  rate  of  the  dissipa¬ 
tion  by  av  =0^/0°,  so  that  the  resulting  production 
rate  Pw  is  basically  unaffected  by  the  WNL  correc¬ 
tion.  This  way,  the  effect  of  the  correction  is  strength¬ 
ened  because  the  production  of  TKE  is  lowered  and 
the  production  of  dissipation  is  not.  This  choice  is 
adopted  here,  but  for  the  compressible  part  of  Pw,  i.e. 
a(Lj/k)Pkc,  which  is  not  divided  by  av  because  it  does 
not  depend  directly  on  /it.  Durbin21  also  recommends 
dividing  the  dissipation  rate  of  the  dissipation  by 
av.  This  is  not  justified  by  any  dependency  of  Dw 
on  nt,  so  that  this  choice  is  not  retained  here,  which 
further  strengthens  the  effect  of  the  WNL  correction. 

Link  with  Non-Linear  Models 

The  link  of  the  Weakly  Non-Linear  correction  (8) 
with  fully  non-linear  models  is  clearly  pointed  out 
in  the  work  by  Pope,28  when  deriving  an  effective 
viscosity  more  general  than  the  linear  Boussinesq  as¬ 
sumption:  the  first  term  in  the  full  development  of  the 
Reynolds  stress  has  the  same  form  as  proposed  here. 

Successfull  Applications  of  WNL  Corrections 

The  preceeding  analysis  shows  that  1)  the  gener¬ 
alization  of  the  empirical  Bradshaw’s  assumption,  2) 
the  enforcement  of  the  realizability  principle,  and  3) 
the  search  for  the  most  general  effective  eddy  vis¬ 
cosity  all  end  up  with  very  similar  dependencies  of 
cM  on  the  strain  rate  and  vorticity  invariants  in  two- 
equation  turbulence  models.  These  expressions  have 
been  successfully  used  to  improve  the  prediction  of 
flows  with  strong  velocity  gradients  such  as  flows  over 
turbine  blades,22  airfoils  in  a  turbulent  freestream21 
and  impinging  jets.24  In  Crossing  SWBLI,  the  pro¬ 
posed  WNL  correction  has  been  shown  to  efficiently 
suppress  the  violation  of  realizability  and  to  improve 
the  prediction  of  wall  pressure  and  heat-transfer  coeffi¬ 
cient,  even  though  further  progress  is  still  needed.17  In 
the  next  section,  the  ability  of  the  correction  to  address 
the  pressure  and  skin-friction  issues  in  3-D  SWBLI  is 
assessed. 

Performance  of  WNL  Models  in  Single 
Fin  SWBLI 

Single  Sharp  Fin  Plate  Experiments 

The  single  fin  plate  interaction  is  representative  of 
flows  developing  over  many  elements  of  supersonic 
bodies,  such  as  inlets,  fuselage/wing  junctions,  con¬ 
trol  surfaces. . .  In  this  configuration  using  a  simple 
geometrical  definition  (Fig.  1),  the  shock  sheet  em¬ 
anating  from  the  fin  leading  edge  interacts  with  the 


incoming  boundary  layer  developing  along  the  bot¬ 
tom  wall.  The  boundary  layer  separates  and  the  foot 
of  the  shock  takes  a  lambda  shape.  Experiments, 
and  then  computations,  have  revealed  the  vortical  na¬ 
ture  of  the  flow  embedded  within  the  lambda  foot. 
The  classification  of  the  secondary  separation  devel¬ 
opment  stages  has  been  established  by  Zheltovodov:25 
types  I  and  II  correspond  to  weak  interactions  without 
secondary  separation;  when  the  interaction  strength 
increases,  a  secondary  separation  first  appears  below 
the  vortical  flow  (type  III),  then  tends  to  disappear 
(types  IV,  V)  and  re-appears  closer  to  the  fin  for 
very  strong  interactions  (type  VI).  The  attachment 
of  the  vortical  flow  on  the  bottom  plate  yields  large 
peaks  in  pressure  and  skin-friction  coefficient.  The 
primary  features  of  the  interactions  (wall  pressure  dis¬ 
tribution,  primary  separation  and  attachment  lines) 
are  fairly  well  predicted  by  RANS  computations.4 
However,  the  skin-friction  peak  is  overestimated,  of¬ 
ten  by  more  than  50  %,  and  no  model  has  been 
yet  shown  to  be  able  to  predict  the  secondary  sepa¬ 
ration  appearance/disappearance/re- appearance  phe¬ 
nomenon.5  Here,  three  configurations  denoted  SF315, 
SF420  and  SF430  (Table  1),  are  selected  to  reproduce 
the  interactions  of  type  III,  IV  and  VI. 

Numerical  Method 

The  three-dimensional  compressible  steady  RANS 
equations  are  solved  using  the  code  GASPex.29  The 
inviscid  fluxes  are  computed  to  third-order  accuracy 
using  the  Roe  scheme  and  a  MUSCL  reconstruction 
method  with  the  Min-Mod  limiter.  Turbulence  is  rep¬ 
resented  using  the  Wilcox’  k—ui  model,30  implemented 
with  the  smooth  wall  boundary  condition31  and  de¬ 
noted  WI.  The  weakly  non-linear  version  of  this  model, 
denoted  WD+,17  is  the  extension  to  compressible  flows 
of  the  non-linear  correction  of  Durbin21  as  described 
above.  The  subsequent  results  are  affected  neither  by 
the  wall  condition  nor  the  freestream  condition  on  the 
k  —  to  model,  in  the  limits  given  in  Ref.  32. 
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Table  1  Single  sharp  fin  plate  flow  features:  Mach  number,  fin  angle,  total  temperature,  stagnation 
pressure,  station  of  the  boundary  layer  measurements,  Mach  number  at  this  station,  wall  temperature, 
unit  Reynolds  number,  boundary  layer  thickness,  momentum  thickness. 


Case 

Moo 

Po 

(deg.) 

T 

-Lioo 

(K) 

p. 

1  ICO 

(kPa) 

A0 

(mm) 

M0 

T 

±  w 

(K) 

Ruo 

(106/m) 

do 

(mm) 

0o 

(mm) 

SF315  (Ref.  25) 

3 

15 

283 

422 

-4 

2.87 

265 

36.56 

5.3 

0.36 

SF420  (Ref.  26) 

3.97 

20 

293 

1524 

38 

3.97 

283 

72.25 

2.87 

0.128 

SF430  (Ref.  27) 

4 

30.6 

283 

610 

0 

3.77 

260 

33.63 

4.9 

0.29 

The  inflow  and  lateral  conditions  correspond  to  a 
previously  computed  flat  plate  boundary  layer  match¬ 
ing  the  momentum  thickness  as  measured  in  the  ex¬ 
periments.  The  computational  domain  is  limited  to 
the  region  influenced  by  the  interaction  (Fig.  2).  The 
longitudinal  step  is  about  one  third  of  the  Incom¬ 
ing  Boundary  Layer  (IBL)  thickness.  In  the  direction 
normal  to  the  bottom  plate,  100  cells  out  of  160  are 
exponentially  distributed  along  the  IBL  height.  The 
first  cell  height  is  0.18  pm  in  the  SF315  case  (resp. 
0.10  pm  and  0.15  pm  in  the  SF420  and  SF430  cases). 
It  represents  less  than  0.05  wall  units  in  the  IBL.  The 
maximum  height  in  wall  units  on  the  bottom  plate  is 
0.24  and  it  is  reached  along  the  attachment  line.  In 
the  transverse  direction,  one  half  of  the  160  cells  is 
concentrated  within  the  side- wall  boundary  layer  with 
a  first  cell  width  of  about  the  double  of  the  first  cell 
height.  The  total  number  of  cells  is  around  3.3  million 
(128  x  160  x  160)  in  the  SF315  and  SF430  cases,  2.8 
million  (112  x  160  x  160)  in  the  SF420  case. 

Solving  the  Grid- Convergence  Issue 

The  mesh-sequencing  technique29  is  used  to  assess 
the  grid-convergence  of  the  solutions.  The  preceeding 
description  corresponds  to  the  very  fine  level  of  the 
grids  ( Nx  x  160  x  160).  The  fine  level  ( Nx  x  80  x  80)  is 
obtained  by  merging  the  cells  two  by  two  in  both  nor¬ 
mal  and  transverse  directions.  Medium  ( Nx  x  40  x  40) 
and  coarse  levels  ( Nx  x  20  x  20)  are  defined  by  repeating 
the  same  procedure.  Steady  solutions  are  obtained  on 
all  four  levels,  from  the  coarser  to  the  finer.  Figure  3 


(a)  P/Pi,  wi  (b)  Cf  x  1000,  wi 


(c)  P/Pi,  wd+  (d)  Cf  x  1000,  wd+ 


Fig.  3  Grid  convergence  for  the  SF430 
configuration  at  the  wall  in  the  cross-section 
Xi  =  122.5  mm. 


shows  the  wall  distribution  of  pressure  (normalized  by 
the  pressure  in  the  IBL)  and  skin-friction  coefficient  in 
cross-section  X\  =  122.5  mm  of  the  SF430  configura¬ 
tion.  It  demonstrates  that  with  the  wi  model,  grid- 
convergence  is  achieved  very  fast  everywhere,  except 
in  the  vicinity  of  the  attachment  line  Z\lX\  ~  0.1.  In 
this  region,  a  constant  growth  in  skin-friction  coeffi¬ 
cient  (Fig.  3b)  and  even  pressure  (Fig.  3a)  is  observed 
when  grid  is  refined.  This  misbehavior  is  cured  by 
the  WNL  correction  as  illustrated  by  Figs.  3c  and  3d 
showing  that  the  solution  on  the  very-fine  grid  is  grid- 
converged  everywhere. 

This  has  to  be  related  to  the  fact  that  some  terms  in 
the  turbulent  equations  of  a  linear  two-equation  model 
tend  toward  infinity  across  the  shock  when  grid  is  re¬ 
fined.  For  instance,  if  is  the  width  of  the  cell 
crossed  by  the  shock  wave,  and  if  Av  is  the  corre¬ 
sponding  finite  jump  in  velocity,  the  gradients  verify: 


S2 


4 

3 


and 


Skk  ~ 


Av 


(12) 


Consequently,  the  asymptotic  behavior  of  the  incom¬ 
pressible  and  compressible  parts  of  the  TKE  produc- 
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(a)  SF315,  wi 


(b)  SF315,  oil-flow  pattern  (c)  SF315,  wd+ 


(d)  SF430,  wi 


(e)  SF430,  oil-flow  visual. 


(f)  SF430,  wd"* 


0  20  40  60  80  100  120  140  160  180  200 

Z,  (mm) 


0  20  40  60  80  100  120  140  160  180  200 

Z,  (mm) 


Fig.  4  Skin- friction  lines  on  the  bottom  plate  of  single  sharp  fin  plate  configurations  (Top  views). 


tion  rate,  integrated  over 

the  cell  volume  is: 

J  PkidV  ~ 

4  c„  Ad 

—  pk  —  ——  Av 

3  oj  A£ 

(13) 

1  PkcdV  ~ 

2 

—  -  pk  Av 

O 

(14) 

With  a  linear  Boussinesq  model  (cM  =  c°),  f  Pk  dV  be¬ 
haves  as  1/A£  in  the  cell  crossed  by  the  shock  wave. 
This  unphysical  behavior  prevents  grid-convergence. 
With  the  WNL  correction  Eq.  (8),  f  Pk  dV  tends  to¬ 
ward  the  finite  value  (2/3)  (c  —  1)  pk  Av  when  grid  is 
refined  (A£  tends  toward  0).  Then,  grid-convergence 
becomes  possible.  In  most  configurations,  the  mis¬ 
behavior  of  the  TKE  production  rate  has  no  practical 
consequence,  because  the  fluid  crossing  the  shock  wave 
does  not  reach  the  regions  under  observation.  How¬ 
ever,  in  the  considered  SWBLI,  the  fluid  originating 
from  the  upper  part  of  the  boundary  layer  crosses  the 
lambda-foot  shock  and  is  entrained  by  the  primary 
vortex  toward  the  bottom  wall  that  it  impinges  along 
the  attachment  line.17  It  becomes  highly  probable  that 
a  flaw  in  grid-convergence  near  the  shock  is  reflected 
in  the  solution  along  the  attachment  line. 


Solving  the  Secondary  Separation  Issue 

Views  of  the  skin-friction  lines  on  the  bottom  plate 
are  reported  in  Fig.  4  in  the  weak  (SF315)  and  strong 
(SF430)  interaction  cases.  The  angles  of  the  attach¬ 
ment  and  separation  lines  with  respect  to  the  incoming 
flow  direction  are  reported  in  Table  2.  The  domain 
under  influence  of  the  interaction  is  bounded  by  the 
primary  separation  line  Si .  Its  angle  is  systematically 
underestimated  by  the  linear  model  by  1.5°  to  4°  and 
overestimated  by  the  wd+  model  by  1°  to  2.5°.  The 
angle  of  the  primary  attachment  line  Ri  is  correctly 


Configuration 

PA 

SS 

SA 

PS 

SF315 

EXP 

23.5 

36.5 

— 

41.5 

WI 

24.0 

— 

— 

40.0 

WD+ 

23.0 

27-31 

— 

42.5 

SF420 

EXP 

25.0 

39.0 

— 

46.0 

WI 

26.0 

— 

— 

44.5 

WD+ 

26.0 

— 

— 

47.5 

SF430 

EXP 

37.0 

46.5 

48.0 

66.0 

WI 

38.8 

— 

— 

62.0 

WD+ 

39.5 

47.0 

50.5 

68.5 

Table  2  Angles  of  Primary  Attachment  (PA), 
Secondary  Separation  (SS),  Secondary  Attachment 
(SA),  and  Primary  Separation  (PS)  lines. 
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(a)  SF315  (b)  SF420  (c)  SF430 


Fig.  5  Wall  pressure  and  skin-friction  coefficient  in  cross-sections  of  single  sharp  fin  plate  configurations. 


predicted  by  both  models.  In  the  type  VI  interaction 
(SF430),  the  secondary  separation  and  re-attachment 
lines  S-2  and  R-2  visible  in  the  experiment  (Fig.  4e)  are 
obviously  predicted  by  the  WD+  model  (Fig.  4f)  while 
they  are  missed  by  the  wi  model  (Fig.  4d).  The  an¬ 
gles  of  these  lines  are  very  close  to  the  measurements 
(Table  2).  Now,  in  the  type  III  interaction  (SF315), 
the  other  kind  of  secondary  separation  as  observed  in 
the  experiment  (Fig.  4b)  is  absent  from  both  numer¬ 
ical  solutions.  Nevertheless,  the  WD+  model  gives  an 
obvious  deviation  of  the  skin-friction  lines  (Fig.  4c), 
which  indicates  a  good  tendency  of  the  model. 

Solving  the  Pressure  and  Skin-Friction  Issues 

The  WNL  correction  has  also  a  beneficial  effect  on 
the  pressure  distribution.  In  Figs.  5a-5c,  the  com¬ 
puted  pressures  on  the  bottom  plate  are  compared 
with  measurements  along  cross-sections  normal  to  the 
fin  surface.  The  cross-section  line  is  straight  in  the 
SF315  and  SF430  cases  (Figs.  5a,  5c)  and  circular  in 
the  SF420  case  (Fig.  5b).  The  geometrical  variables 
Xi,  Z\ .  R  and  B  are  defined  in  Fig.  1.  The  predic¬ 
tion  of  the  pressure  plateau  under  the  vortical  flow 
embedded  within  the  lambda  foot  of  the  shock  is  poor 
without  the  WNL  correction  and  very  good  with  it 


(SF315:  Zi/A'i  =  0.35  -  0.55,  SF420:  B  =  35  -  50°, 
SF430:  Z\jX\  =0.35  —  0.8).  Figures  5a-5c  confirm 
that  the  interaction  extent  (the  region  where  the  pres¬ 
sure  differs  from  its  level  in  the  IBL)  is  overestimated 
with  the  WNL  model  and  underestimated  without  it. 

The  local  dip  in  pressure  observed  between  the  at¬ 
tachment  and  the  vortical  flow  is  well  predicted  by 
the  WD+  model  while  it  is  missed  by  the  linear  model 
(SF315,  SF430:  Zl/Xl  ~  0.3,  SF420:  B  -  32°).  In 
the  strong  interaction  case,  the  WD+  model  gives  the 
secondary  shock  and  the  weak  local  peak  in  pressure 
associated  with  the  secondary  re-attachment  (Fig.  5c, 
Zi/Xi  ~  0.35). 

Finally,  the  pressure  distribution  between  the  fin 
and  the  primary  attachment  line  is  nearly  the  same 
in  both  numerical  solutions.  The  peak  value  of  the 
pressure  is  underestimated  by  about  9%  in  the  mild 
interaction  (SF315),  and  overestimated  by  about  5% 
in  the  strong  interactions  (SF420  and  SF430) . 

Figures  5d-5f  show  the  skin-friction  coefficient  Cf  in 
the  same  cross-sections  as  for  the  wall  pressure.  The 
major  effect  of  the  WNL  correction  is  to  decrease  the 
maximum  in  Cf  by  30  to  35  %.  In  the  SF420  case, 
two  series  of  measurements  are  reported  (Fig.  5e).  The 
second  one  (closed  circles  labelled  EXP-2)  resulted  in 
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corrected  values  of  the  skin-friction  coefficient  near  the 
primary  attachment  line.4  The  WD+  solution  agrees 
well  with  the  correct  measurements.  The  wi  model, 
like  other  linear  two-equation  models4  or  even  the 
Baldwin  and  Lomax  model  as  modified  by  Panaras,6 
yields  an  overestimate  by  at  least  50%  of  the  maximum 
skin-friction  coefficient.  Some  other  models  of  k  —  e 
type  do  not  exhibit  this  erroneous  peak  value,  but  the 
overall  Cf  distribution  is  quite  far  from  the  measure¬ 
ments  (Ref.  4,  Fig.  50).  Although  measurements  are 
not  available  in  the  SF430  case,  the  reduction  in  Cf 
near  Z\jX\  ~  0.35  (Fig.  5f)  is  validated  by  the  sub¬ 
sequent  appearance  of  the  secondary  separation  and 
re-attachment  lines  (Fig.  4e,  4f). 

Reducing  the  Turbulent  Kinetic  Energy  Level 

Figure  6a  shows  the  WNL  correction  factor 
av  =  cM/c°  'n  cross-section  X\  =  122.5  mm  of  the 
SF430  configuration.  The  WNL  correction  is  active 
(i.e.  av  is  smaller  than  1)  only  through  the  legs  of  the 
shock  wave  and  also  in  the  core  of  the  vortical  flow  and 
near  the  attachment  line.  The  effect  on  the  turbulence 
intensity  k/U is  dramatic  :  the  maximum,  achieved 
in  the  core  of  the  vortex,  is  more  than  halved  by  the 
WNL  correction  (Figs.  6b,  6c).  Consequently,  the  in¬ 
tensity  of  the  rear-leg  shock  is  larger,  the  secondary 

(a)  wd+,  WNL  factor  av  =  cM/c° 


(b)  WI,  turbulence  intensity  k/U ^ 


(c)  WD+,  turbulence  intensity  k/U ^ 


Fig.  6  Turbulent  quantities  in  cross-section 
Xi  =  122.5  mm  of  SF430  configuration. 


shock  can  appear  and  the  primary  vortex  is  flatter. 
The  action  of  the  WNL  correction  through  the  fore¬ 
leg  of  the  shock  yields  the  observed  wider  extent  of 
the  interaction.  Due  to  the  sensitivity  of  the  solution 
in  this  region,  only  a  very  slight  change  in  the  correc¬ 
tion  should  be  sufficient  to  improve  the  prediction  of 
the  interaction  extent. 

Performance  of  WNL  Models  in 
Transonic  Flows 

The  identified  WNL  correction  has  been  successfully 
used  in  supersonic  SWBLI.  It  was  initially  designed 
to  handle  subsonic  configurations  with  highly  curved 
streamlines  around  stagnation  points.  In  order  to 
complete  the  validation,  transonic  cases  have  to  be 
considered  as  well.  The  purpose  of  this  section  is  to 
assess  the  behavior  of  the  WNL  correction  in  two  dif¬ 
ferent  transonic  SWBLI.  The  second  configuration  is 
chosen  to  give  a  flavor  of  the  interest  of  weakly  non¬ 
linear  models  for  unsteady  flows. 

Transonic  SWBLI  in  a  Channel  Flow 

The  first  transonic  configuration  is  a  two- 
dimensional  channel  flow  with  a  bump  on  the  lower 
wall  (Fig.  7).  The  channel  inlet  Mach  number  is  0.615 
with  a  total  pressure  of  96  kPa  and  a  total  temperature 
of  300  K;  the  back  pressure  is  set  to  61.5  kPa.  The  ex¬ 
perimental  study  due  to  Delery33  shows  that  along  the 
bump  ramp,  the  flow  accelerates  and  reaches  a  Mach 
number  of  about  1.4.  A  shock  forms  to  accomodate 
the  back  pressure.  On  the  lower  wall,  above  the  down¬ 
stream  part  of  the  bump,  the  shock- wave/boundary- 
layer  interaction  causes  the  separation  of  the  bound¬ 
ary  layer,  the  formation  of  a  recirculation  bubble  and 
a  large  lambda-foot  shock  (Fig.  8a).  This  configu¬ 
ration  has  been  considered  in  several  computational 
workshops,35,36  from  which  it  results  that  no  linear 
two-equation  model  nor  tested  Reynolds-Stress  model 
is  able  to  reproduce  the  correct  large  lambda  shock 
visible  in  the  experiment  nor  the  pressure  plateau  un¬ 
derneath.  Only  the  SST  k  —  lj  model  and  the  cubic 
Non-Linear  Eddy-Viscosity  Model  (NLEVM)  by  Craft 
et  al.37  show  some  ability  to  exhibit  these  features.36 


Top  wall 

r=20 

\ 

)-0.J5  * 

Z 

Circular  r  =  422.67 

'O 

Ov 

r=20 

V 

O 

t 

156.88 

29.49 

100 

X 

286.37 

Dimensions  in  mm 

Fig.  7  Geometry  of  the  transonic  bump33  (from 
Ref.  34,  with  permission). 
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In  a  recent  numerical  study  discussed  hereafter, 
Dieudonne34  contributes  to  a  better  understanding  of 
the  situation.  Across  both  legs  of  the  lambda  shock, 
the  TKE  production  term  given  by  the  usual  linear 
stress-strain  relation  grows  unbounded.  Consequently, 
the  eddy  viscosity  is  too  high  and  prevents  the  devel¬ 
opment  of  the  lambda  shape  of  the  shock.  Figure  8b 
shows  a  typical  result  from  linear  two-equation  mod¬ 
els.  This  one  is  obtained  with  the  wi  k  —  u)  model. 
By  adding  the  SST  correction,  a  much  better  lambda 
shape  is  obtained  (Fig.  8c).  Dieudonne  further  stud¬ 
ied  this  case  by  considering  first  the  cubic  NLEVM  of 
Craft  et  al.37  coupled  with  the  k— lj  model.  The  results 
are  even  worse  than  with  the  linear  model.  Nonethe¬ 
less,  by  rising  the  back  pressure  by  5%  and  reducing 
the  eddy  viscosity  by  3/4  in  the  cubic  NLEVM,  a  very 
good  pressure  distribution  and  lambda  shape  are  ob¬ 
tained  (Fig.  8d),  with  a  slight  improvement  over  the 
SST  results.  Second,  Dieudonne  tested  the  current 
WD+  model  with  no  adjustment  and  reproduced  the 
better  results  as  obtained  with  the  ad-hoc  NLEVM: 
the  lambda  shape  of  the  shock  (Fig.  8e)  is  comparable 
with  the  one  given  by  the  ad-hoc  NLEVM  (Fig.  8d), 

(a)  Experiment  (b)  WI 


(c)  SST 


(d)  ad-hoc  cubic  NLEVM 


(e)  wd+ 


(f)  Exp.  /  WD+ 


Fig.  8  Transonic  bump  flow  in  the  separation  re¬ 
gion;  (a):  experimental  interferogram,  (b— e):  com¬ 
puted  isobars,  (f):  density  isolines  (wd4  model) 
with  superimposed  lines  from  the  experimental  in¬ 
terferogram  (from  Ref.  34,  with  permission). 


and  a  superimposition  of  the  density  field  with  lines  ex¬ 
tracted  from  the  experimental  interferogram  (Fig.  8f) 
shows  a  very  good  agreement.  Further  details  about 
these  results,  including  pressure  distributions  and  ve¬ 
locity  profiles  can  be  found  in  Ref  34. 

This  work  supports  the  idea  that,  in  the  case  of  the 
transonic  bump,  the  key-point  is  to  lower  the  eddy  vis¬ 
cosity,  and  consequently  the  TKE  level  in  the  separat¬ 
ing  flow,  rather  than  accurately  modeling  the  Reynolds 
stress  anisotropy.  The  latter  way  could  certainly  fur¬ 
ther  improve  the  results,  but  much  of  the  work  is  done 
by  the  WNL  correction. 

Unsteady  Transonic  SWBLI  over  an  Airfoil 

As  underlined  by  Dolling,1  RANS  computations 
with  standard  linear  two-equation  models  are  known 
to  be  unable  to  predict  pressure  fluctuations  when  the 
SWBLI  causes  the  appearance  of  unsteadiness  in  the 
flow.  An  archetype  of  this  phenomenon,  which  is  of 
primary  importance  for  aircraft  wing  design,  is  the  ap¬ 
pearance  of  flow  oscillations  induced  by  SWBLI  over 
a  transonic  airfoil.  Depending  on  the  airfoil  shape 
and  flow  conditions,  Shock-Induced  Oscillations  (SIO) 
of  the  flow  may  appear,  coupled  with  periodic  exten¬ 
sion/reduction  of  the  separated  region  downstream  the 
shock,  and  they  can  yield  airfoil  buffeting  through  the 
mechanical  response  of  the  wing  structure. 

Several  experimental  studies  have  been  carried  out 
to  understand  the  physics  of  this  phenomenon.  McDe- 
vitt  &  Okuno38  have  experimentally  identified  the 
Mach-number /incidence- angle  relationship  for  the  SIO 
onset  over  the  NACA0012  airfoil  at  Mach  numbers 
between  0.7-0. 8,  incidence  angles  between  0-5°  and 
chord  Reynolds  number  between  1-14  106.  Barakos  & 
Drikakis39  have  shown  that  no  combination  of  Mach 
number  and  incidence  angle  in  this  range  can  yield  SIO 
using  standard  linear  k  —  e  turbulence  models.  On  the 
opposite,  either  the  k  —  e  or  the  k  —  ui  version  of  the 
cubic  NLEVM  by  Craft  et  al.37  does  a  good  job  in 
predicting  the  Mach-number /incidence-angle  relation¬ 
ship  for  the  SIO  onset.  Barakos  &  Drikakis  have  also 
established  that  by  suppressing  the  functional  depen¬ 
dence  of  cM  on  the  strain  and  vorticity  invariants  in  the 
NLEVM  model,  no  more  SIO  can  be  obtained.  They 
conclude  that  the  key-point  in  the  success  of  both  cu¬ 
bic  NLEVM  is  the  functional  cM  and  not  the  non-linear 
expansion  of  the  shear  stress,39  so  that  it  would  be 
worthwhile  to  use  a  functional  in  conjunction  with 
a  linear  two-equation  model. 

On  the  other  hand,  in  the  ONERA  S3MA  wind 
tunnel,  unsteady  pressure  and  skin-friction  coefficients 
have  been  measured  on  the  RA16SC1  supercritical 
airfoil  in  a  flow  at  a  Mach  number  of  0.732  for  a 
chord  Reynolds  number  of  4.2  106  at  different  angles 
of  attack  a,,  .40  SIO  have  been  measured  only  for  ae 
between  3-5°.  Inviscid/ viscous  coupled  computations 
indicate  that  the  effects  of  the  side- wall  boundary  lay- 
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ers  on  the  symmetry  plane  may  be  taken  into  account 
by  performing  the  two-dimensional  computations  at 
lower  Mach  number  (0.723)  and  angle  of  attack  (ac  = 
ae  —  1°).40  In  the  case  ae  =  4°  corresponding  to  the 
maximum  amplitude  in  the  flow  oscillations  (Fig.  9) 
and  in  the  measured  pressure  fluctuations,  Furlano  et 
al .41  show  that  no  oscillation  can  be  obtained  using 
the  Smith  k  —  t  turbulence  model42  and  that  pressure 
oscillations  appear  with  correct  frequency  and  ampli¬ 
tudes  when  the  k  —  l  linear  eddy- viscosity  is  arbitrarily 
limited  to  an  ad-hoc  maximum  value.  A  similar  phe¬ 
nomenon  is  observed  with  a  k— e  model43  by  Rouzeaud 
et  al.AA  These  results  support  the  idea  that  the  key- 
point  is  the  limitation  of  the  eddy-viscosity  level  as 
obtained  with  a  linear  turbulence  model. 

A  recent  numerical  study  contributes  to  assert  this 
hypothesis.  Goncalves  &  Houdeville45  have  assessed 
the  effect  of  a  WNL  correction  coupled  either  with 
the  high-Reynolds  number  k  —  e  or  with  the  k  —  lj 
turbulence  models.  Wall  effects  are  modeled  through 
a  wall-law  approach.46  In  the  considered  case  of  the 
RA16SC1  airfoil  at  ae  =  4°  as  described  above,  it 
has  been  verified  that  the  choice  of  the  SST  correc¬ 
tion  or  the  WNL  correction  Eq.  (8),  both  coupled 
with  the  similar  transport  equations  (9-10),  do  not 
really  matter.  This  is  understandable  because  in  a 
two-dimensional  boundary  layer,  both  invariants  s  and 
w  degenerate  in  the  same  limit. 

Goncalves  &  Houdeville45  show  that  with  the  linear 
turbulence  models  (i.e.  without  the  SST  correction), 
no  unsteady  fluctuations  are  obtained  with  the  k  —  lu 
model,  and  some  develop  with  the  k  —  e  model.  In  the 
latter  case  however,  Fig.  10a  shows  that  the  exten¬ 
sion  of  the  fluctuating  region  and  the  amplitude  of  the 
pressure  fluctuations  are  clearly  underestimated.  The 
same  figure  shows  the  dramatic  and  beneficial  effect  of 
the  SST  correction  with  both  basic  turbulence  models: 
the  extension  of  the  fluctuating  region  and  the  ampli¬ 
tude  of  the  pressure  fluctuations  are  fairly  predicted 
by  both  k  —  e  and  k  —  u>  models  with  the  WNL  correc¬ 
tion.  Recent  computations  with  a  low- Reynolds  num¬ 
ber  k  —  e  model43  and  a  dual-time  stepping  method44 
confirm  that  the  conclusion  does  not  depend  on  the 
wall-law  approach.47 

Between  ae  =3  —  5°,  the  SST  k  —  lj  model  predicts 
also  the  measured  evolution  of  the  frequency  of  the 
oscillations  while  this  evolution  is  clearly  missed  with 
the  k  —  e  model  (Fig.  10b).  Finally,  only  the  SST  fe-cu 
model  is  able  to  predict  the  observed  damping  of  the 
oscillations  for  ae  =  5°  (Fig.  10c). 

Conclusions 

Usual  two-equation  turbulence  models  fail  in  pre¬ 
dicting  SWBLI  mainly  because  of  the  linear  relation¬ 
ship  between  the  turbulent  stress  and  the  strain  tensor 
through  the  constant  value  of  in  the  definition  of  the 
eddy  viscosity.  Weakly  Non-Linear  turbulence  models 


Fig.  9  RA16SC1  airfoil,  SST  k  —  oj  model,  Mach 
isolines  at  two  instants  in  phase  opposition  (from 
Ref.  45,  with  permission). 


Fig.  10  RA16SC1  airfoil,  (a):  4°  angle  of  attack, 
root-mean-square  pressure  fluctuations  upstream 
and  on  the  suction  side  of  the  airfoil,  (b):  Strouhal 
number  vs  angle  of  attack,  (c):  amplitude  of  the 
lift-coefficient  fluctuations  vs  angle  of  attack  (from 
Ref.  45,  with  permission). 

can  be  obtained  very  simply  by  considering  as  a 
functional  of  the  strain  rate  and  vorticity  invariants. 
At  least  three  different  theoretical  approaches  all  end 
up  with  a  similar  expression  of  cM:  the  generalization 
of  the  empirical  Bradshaw’s  assumption,19’20  the  en- 
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forcement  of  the  realizability  principle,21’22  and  the 
search  for  a  more  general  effective  eddy  viscosity.28 

Such  a  Weakly  Non-Linear  correction  has  been  pro¬ 
posed  and  applied  to  three  single  sharp  fin  plate  con¬ 
figurations  with  different  SWBLI  strengths.  The  cor¬ 
rection  yields  a  much  better  prediction  of  the  pressure 
under  the  vortical  flow  formed  within  the  lambda  foot 
of  the  shock.  The  maximum  of  the  skin-friction  coeffi¬ 
cient  is  decreased  by  one  third,  yielding  levels  compa¬ 
rable  to  the  measurements.  In  the  strong  interaction 
case  (Mach  4,  30.6°  fin),  the  secondary  separation  due 
to  the  formation  of  a  secondary  shock  appears,  while  it 
is  totally  absent  from  the  solution  without  the  correc¬ 
tion.  The  turbulence  intensity  is  more  than  halved  in 
the  vortical  flow,  which  explains  the  decrease  in  skin- 
friction.  It  is  experienced  that  when  grid  is  refined, 
the  levels  of  pressure  and  skin-friction  coefficient  do 
not  converge  at  the  attachment  with  the  linear  model, 
and  do  converge  with  the  WNL  correction.  It  is  re¬ 
lated  to  the  grid-divergence  of  some  terms  in  the  linear 
model  equations,  which  is  cured  by  the  correction. 

Unlike  the  SST  correction,20  the  WNL  correction 
ensures  realizability  and  is  active  in  a  pure  irrotational 
strain  (like  when  crossing  a  normal  shock  wave  for 
instance).  Nonetheless,  the  two-dimensional  limit  of 
the  present  WNL  and  the  SST  corrections  are  nearly 
the  same.  In  a  two-dimensional  transonic  channel  flow 
over  a  bump,  the  WNL  correction  predicts  the  lambda- 
shape  of  the  shock  at  least  as  well  as  the  SST  correction 
and  an  ad-hoc  version34  of  a  cubic  NLEVM.37  In 
transonic  flows  over  airfoils,  observed  Shock-Induced 
Oscillations  cannot  be  properly  obtained  with  linear 
turbulence  models,  as  shown  by  Barakos  &  Drikakis39 
with  k  —  e,  Furlano  et  al ,41  with  k  —  l,  Goncalves  & 
Houdeville45  with  k  —  u>.  In  all  three  cases,  SIO  can 
be  obtained  either  with  the  cubic  non-linear37  version 
of  the  k  —  e  model,39  an  arbitrary  maximum  limit  on 
the  eddy  viscosity41  or  the  use  of  the  WNL  or  SST 
correction.45  Moreover,  the  WNL  or  SST  k  —  u>  model 
is  able  to  predict  the  evolution  of  the  frequency  and 
amplitude  of  the  oscillations  when  varying  the  angle  of 
attack  of  the  airfoil.45  These  results  confirm  the  state¬ 
ment  by  Barakos  &  Drikakis:39  in  order  to  predict  SIO 
over  airfoils,  the  key-point  is  the  functional  and  not 
the  non-linear  expansion  of  the  shear  stress. 

The  overall  conclusion  is  that  the  WNL  correc¬ 
tion,  originally  developped  for  two-dimensional  sub¬ 
sonic  flows21’22’24  has  been  successfully  used  in  3-D 
supersonic  and  2-D  transonic  steady  and  unsteady 
flows.  Some  flaws  still  persist  in  the  mean  flow  and 
should  be  further  reduced  by  refining  the  functional 
for  intermediate  values  of  the  strain  and/or  vorticity 
invariants,  and  by  introducing  either  a  variable  turbu¬ 
lent  Prandtl  number  or  specific  equations  to  further 
deal  with  the  heat-transfer  issue.  Then,  non-linear 
turbulence  models  could  be  considered  for  a  better  pre¬ 
diction  of  turbulent  stresses. 
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